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ABSTRACT 

A new Monte-Carlo algorithm for calculating time-dependent radiative-transfer 
under the assumption of LTE is presented. Unlike flux-limited diffusion the method 
is polychromatic, includes scattering, and is able to treat the optically thick and free- 
streaming regimes simultaneously. The algorithm is tested on a variety of 1-d and 
2-d problems, and good agreement with benchmark solutions is found. The method 
is used to calculate the time- varying spectral energy distribution from a circumstellar 
disc illuminated by a protostar whose accretion luminosity is varying. It is shown that 
the time lag between the optical variability and the infrared variability results from a 
combination of the photon travel time and the thermal response in the disc, and that 
the lag is an approximately linear function of wavelength. 

Key words: Radiative transfer - methods: numerical - stars: pre-main-sequence - 
planetary systems: protoplanetary discs 



1 INTRODUCTION 

Theoretical astrophysics often necessitates computing the 
energy balance between gas and an ambient radiation field. 
Most frequently this involves solving the equation of radia- 
tive equilibrium, and a wide variety of theoretical tools have 
been developed to tackle this problem efiiciently, ranging 
from those which make simplifying assumptions both to the 
geometry (spherical or plane-parallel) and to the underlying 
microphysics (the gray approximation for example) through 
to a full poly-chromatic treatment in multiple dimensions. 

Although most radiation transfer (RT) codes assume 
thermal balance many interesting phenomena occur out 
of equilibrium, and here a time-dependent approach is re- 
quired. A spectacular example of this occurs during su- 
pernova explosions, and progress is being made in time- 



(20091; 



dependent modelling of such objects e.g. Jack et al, 
[Kromer fc STm|p009t . 

Time-dependent methods are also required for radia- 
tion hydrodynamics (RHD). Broadly speaking the trans- 
port methods can be split into those concerned with ion- 
izing radiation, which generally adopt ray-casting methods 
(e.g. [Dale et al.|2007[[Mac Low et al.|2007| l, and those that 
deal with transport in dusty media which use either short- 
characteristic methods, Monte-Carlo techniques, or the dif- 
fusion approximation (e.g. Hofner et al.|2063 Woitke|2006 
Turner & Stone 20011. RHD codes may be further sub- 



divided into those which are applicable to situations in 
which the radiation transport /thermal equilibrium timescale 



for the gas is shorter than characteristic hydrodynamical 
timescale, and those in which the hydrodynamical timescale 
is the shorter. In the former case the radiation transport may 
be conducted as a sequence of pseudo-equilibrium steps (e.g. 



Hofner et al.|2003||Woitko 2006; Dale et al.|2007||Mac Low| 


et al. 


2007 


Freytag & H6fner||2008 


Acreman et al. 


2010 


1. 



However when the thermal timescale becomes comparable 
to, or indeed larger than, the hydrodynamical timescale 
the time-dependent form of the radiation-transfer equations 
must be solved (e.g. [Krumholz et al.||2007| ). 

In many situations RT is implemented assuming energy 
transport occurs via radiative diffusion of photons. This is a 
good assumption for optically thick regions, but breaks down 
in the optically thin regime, since the mean-free-path of pho- 
tons (and thus the speed of the diffusing radiation field) can 
become arbitrarily large; a fiux-limiter is adopted to control 
this problem. Flux-limited diffusion (FLD) is normally im- 



plemented in the gray approximation ( 
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1, although comparisons with a polychromatic 



treatment demonstrate that this is often quite a poor ap- 
proximation ( jPreibisch et al.||1995[ ). Opaque obstacles also 
cause problems for FLD, since the radiation field can diffuse 
around the obstacle when in fact a shadow should be cast. 

A modified version of the Monte Carlo (MC) radiative 
equilibrium method developed by Lucy ( 1999 1 offers an at- 



tractive route to full time-dependent RT, since it would al- 
low a full polychromatic treatment within which multiple 
scattering may be incorporated. The method is computa- 
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tionally demanding, but it has the advantage of being ex- 
tremely efficient to paraUehze. Indeed the feasibifity of using 
such a method has aheady been demonstrated, at least in 



the case where the thermal timescales are short (|Acreman 
et aL|2010| ) 



Here I present an algorithm for calculating the matter 
and radiation field energy densities, as a function of time, 
for an arbitrary distribution photon of sources embedded in 
an arbitrary distribution of gas under LTE conditions. I de- 
scribe one-dimensional benchmark tests of the method, and 
detail its incorporation into the TORUS radiation-transfer 
code. I test the new version of TORUS against a two- 
dimensional disc benchmark ( Pascucci et al.|[2004 |. Finally 
I give a simple application of the code demonstrating the 
effect of a varying accretion luminosity from a protostar on 
its surrounding dusty circumstellar disc. 



2 METHOD 

We describe the physical quantities (densities, temperatures 
etc) on a cell-centred grid. For a gas in LTE at temperature 
T the rate at which it emits energy is given by 



E — Atv ki,Bu dv, 
'o 



(1) 



where fc^ is the absorption coefficient and is the Planck 
function. The rate at which the same gas absorbs energy is 
given by 



A — 4-jY kvJi, dv, 

10 



(2) 



where is the mean intensity of the radiation field. Clearly 
if the gas is in radiative equilibrium then A — E and we find 

T=(^y (3) 



y AaKp 

where Kp is the Planck-mean absorption coefficient. How- 
ever if we consider gas that is not in radiative equilibrium 
then the nett change in energy density of the gas 



iig = A - E 



(4) 



and, mutatis mutandis, the rate of change in energy density 
of the radiation field is 



iir ^ E - A 



(5) 



Now we consider a gas of volume V at time t. Within this 
volume is a star of luminosity L* . The luminosity of the gas 
is given by 



(6) 



We assume that the temperature of the gas is constant over 
a single timestep At. During this timestep we assume that 
the gas and the star produce Ng and A'^, new photon packets 
respectively. The individual photon packet energies are given 
by 

L„At LtAt 



N, 



(7) 



The energy density is related to the temperature by 



RTp 

(7 - 1)m 



(8) 



where 7? is the gas constant, p is the mass density, 7 is the 
ratio of specific heats and /i is the mean molecular weight. 
We follow |Lucy| ( |1999[ ) and use the result that the energy 
density of the radiation field in the interval {u, v + dv) is 
given by 



■.y = A-KJvdv/c 



(9) 



A photon packet moving between events (scatterings, ab- 
sorptions, crossing grid-cell boundaries) contributes an en- 
ergy for a time l/c (where I is the distance between 
events) to the local energy density. The photon energy den- 
sity is therefore 



1 1 1 



(10) 



where the summation is over all photon packets. Now com- 
bining equations [2] and [9] with equation [lO] we obtain an 
expression for the energy absorption rate: 



A = 



(11) 

The new energy density of the gas may then be calculated 



Ug- ^Ug +(A~ E)At 



This explicit integration scheme will require a careful choice 
of the timestep At, which must be short enough to ensure 
stability while also being long enough that the computation 
remains tractable. Timescale considerations are discussed in 
section 12. II 

It is worth noting that equation [l2] is not the only 
method for updating the internal energy of the gas. Since 
we follow the radiation field in detail we can, for each cell, 
calculate the energy in photon packets entering and leaving 
the cell, and hence the change in radiation energy density 
due to the RT (Ait^^trans)- MC estimators exist for Ur for 
both the start and the end of the timestep, so using 



n + l 



= «" + (£- A) At + AUr 



(13) 



12 



we can find [E — A) At and substitute it into equation 
bypassing the need to find the an estimator for A. In prac- 
tise we found that this method required a larger number of 
photon packets (in order to improve the estimate of Ur and 
Attr.trans) than that needed using an estimator for A. 

A single timestep encompasses of a loop over photon 
packets, each with an individual energy and frequency u. 
Information on photon packets that are 'in flight' at the end 
of a timestep are stored on a stack, to be processed as part 
of the RT during the subsequent timestep. The total number 
of photon packets is the sum of the number of packets on the 
stack {Ns), and the number of new photon packets generated 
during the timestep A'p. The path of each photon packet is 
followed as it propagates through the gas, during which time 
the packet may be scattered multiple times. The random 
walk of the photon packet ceases when it is either absorbed 
or leaves the computation domain (in which case the photon 
packet is destroyed), or when its flight time il/c) reaches At 
(in which case the photon packet is added to the stack). Once 
the random walks of all the photons have been calculated, 
we may then use our MC estimate of the absorption rate 



(equation 11 1 to update the gas energy density of each cell 
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via equation |12[ One can immediately see that when the 
photon flight time of a mean free path Kp/c ^ At then 
Na will dominate over A'j, as the calculation proceeds, but 
eventually if radiative equilibrium is reached there will be an 
approximately constant number of packets at each timestep. 



It should be noted that, unlike the Lucy ( 1999 \ radiative 



equilibrium algorithm, energy is not implicitly conserved 
here, as photon packets are created and destroyed during 
each timestep. The quality of the energy conservation is con- 
trolled by the accuracy of our estimator of the absorption 
rate (which in turn is dependent on our MC statistics), and 
by our choice of At. Energy conservation is addressed in 
section |3l 

The algorithm itself is described as a sequence of steps: 

(i) The current energy density of the gas (ug) is used to 
compute the temperature distribution (equation|8| and thus 
the gas energy emission rate as given by equation nj A prob- 
ability distribution of a photon packet being emitted is cal- 
culated over all cells 



Pi 



(14) 



where Vi is the volume of the ith cell. 

(ii) The probability of a photon being produced in the 
gas is calculated by 

La 



Pa 



(15) 



We take r) to be the uniform random deviate. If 77 < pg then 
the photon packet is emitted within the gas, and its position 
is found according to the probability distribution given by 
equation 1 14| otherwise the photon is produced by the star. 

(iii) If a photon is produced by the gas then the photon 
frequency is assigned according to the equation 

V ^ ji^dv/ I j^dv (16) 
Jo Jo 

where jV = kvB,,. If the photon packet is stellar in origin 
then its frequency comes from 



77 = 



B^dv/ 



(17) 



where Bi, is the Planckian at the stellar effective temper- 
ature (although other stellar spectral energy distributions 
may be easily incorporated). 

(iv) A photon can propagate an optical depth 

r = (fc, + fcsca,.)p^ (18) 

(where fcsca.i/ is the scattering coefficient at frequency i^) 
before it is either absorbed or scattered. The optical depth 
is selected randomly: 

r = -ln(l-?7) (19) 

If r is greater than the optical depth to the cell boundary 
then the photon packet position is moved to the cell bound- 
ary and a new optical depth is selected. If the time to travel 
optical depth r is greater than At then then photon packet is 
moved by distance cAt and the photon information is placed 
on the stack. Otherwise the photon position is updated and 
the type of interaction is determined by the albedo 



If the photon packet is scattered then the new direction vec- 
tor is chosen randomly from the appropriate phase function. 
If the photon is absorbed then the appropriate path length 
is stored and the next photon packet is selected. 

(v) Once all photons packets have been processed the new 
absorption rates are computed and the energy density of the 
gas is updated using equation [4] Subsequently the duration 
of the next timestep is determined (see below). 

2.1 Timescales 

Cells in which the emission rate exceeds the absorption rate 
are cooling, and the timescale for cooling is calculated as 



^cool — 



E-A 



(21) 



Conversely, if the energy absorption rate exceeds the emis- 
sion rate for a given cell then its temperature is increasing. 
We may then calculate an approximate radiative equilibrium 
timescale by first computing the equilibrium energy density 
Me (using equation [3| and then finding 



tc 



A~ E 



(22) 



Cells for which tcq < At have their temperatures set to the 
radiative equilibrium temperature. Note that this formalism 
assumes that A and E are constant over AT, and care must 
be taken to ensure that one does not erroneously set cells to 
their equilibrium temperature due to a poorly chosen AT. 



3 TESTS 

I have implemented a number of test cases in order to ex- 
amine the efficacy of the new algorithm. 



3.1 Radiative equilibrium 

If we immerse an absorbing gas in a radiation field of much 
higher energy density the gas will eventually come into equi- 
librium with the radiation. Since Ug <^ Ur we can assume 
Ur is a constant and the evolution of Ug can then be found 



from the ordinary differential equation 
dUn 



dt 



A'KKB{Ug). 



(23) 



For this t est w e adopted the parameters given by |Turner| 



& Stone 



(2001|: p = 10"' gem 
5/3 and Ur = 



lO^'^ erg cm 



= 4 X 10"° cm"\ 
^. We ran two 



p = 0.6, 7 

separate cases, one with Ug intially well below the equilib- 
rium value {ug = 10^ ergcm"^) and one with it well above 
{ug — 10^° ergcm""^). The results are plotted in Figure [l] 
and the agreement with the analytical solution is excellent. 

The second test involved the same physical parameters 
for the box detailed above, but instead of filling the box 
with a high density radiation field we instead set Ur ~ 
and Ug — 10* everywhere. As the RT is followed the mate- 
rial in the box will cool and emit radiation. Since we have 
refiective boundary conditions for the photon packets the 
gas and radiation field will eventually settle into radiative 
equilibrium. We set Np to 1, which means that the gas emits 
one photon per timestep. In Figure[2]we plot the evolution of 
the radiation and gas energy densities with time, and these 



4 Tim J. Harries 




1x10-'^ 1x10-'^ 1x10-" 1x10-' 1x10-' 

Time (seconds) 



Figure 1. The results of the first radiative equilibrium test case 
described in Scction |3.1| In the upper panel the Monte-Carlo re- 
sults are shown as dots, while the analytical solution found from 
equation |23| is displayed as a solid line. The bottom panel shows 
the fractional difference between the analytical and numerical so- 
lutions for the energy density. 
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sho-w good agreement with a simple numerical integration 
of equations |5] and [4] We find that the total energy is con- 
served to better than 0.5% throughout the duration of the 
MC calculation. 

The above tests suggests that the physics of matter- 
radiation interaction is adequately captured by the new al- 
gorithm, and we now progress to tests of the transport of 
radiation. 



3.2 Diffusion limit 

In the optically-thick limit the radiation transport occurs 
according to the diffusion equation 

j2„ 



dUr 

~dt 



-D 



d Ur 

dx^ 



where D is the diffusion coefficient, given by 



(24) 



(25) 



Here we test the Monte-Carlo algorithm on a heat kernel 
problem, in which energy is deposited at a point at t = s 
and is allowed to diffuse through the medium. The analyt- 
ical solution to the diffusion equation for this scenario is a 
gaussian 



u{x, t) — 



exp 



ADt I 



(26) 



We assume a f cm one-dimensional box, divided into lOf 
evenly-spaced bins, with reflective boundary conditions for 
the photons (i.e. adiabatic). We adopt k = fO^'^cm^^ and 
deposit fO^° of energy in photons into the central bin at t = 
Os. (The radiation field immediately comes into radiative 
equilibrium with the material in the box, but Ur 3> ite). 
We follow the radiation transport using the Monte-Carlo 
algorithm, and also by solving equation [24] using the Crank- 
Nicolson method. 

We initially assume pure scattering [a = 1), in which 
case the MC algorithm will always perfectly preserve en- 



Figure 2. The results of the second radiative equilibrium test 
case described in Section [3.1| In the upper panel the Monte-Carlo 
results are shown as solid lines, while the solutions found from 
numerical integration of equations [5] and |4] are displayed as dots. 
The middle and lower panels shows the fractional difference be- 
tween the energy densities from the analytical and Monte Carlo 
solutions for equations [s] and [4] respectively. 



ergy since no photons are created or destroyed and mat- 
ter/radiation energy transfer terms are by definition zero. 
This scenario tests the random walk and photon flight time 
limit section of the MC algorithm, and we flnd good agree- 
ment between the diffusion approximation and the MC al- 
gorithm for this test case (Figure [3|. 

The purely absorptive case (a = 0) is a much more 
challenging proposition. Both the emissivity and absorption 
rates are very high, and Monte-Carlo sampling errors can po- 
tentially lead to deviations from energy conservation. How- 
ever we find good agreement (Figure |4| with the diffusion 
approximation, and energy is conserved to within 2 percent 
even after more than 4000 time steps. 



3.3 Time varying source 

It appears that the diffusion limit is well-modelled using 
the MC algorithm, but such conditions are naturally more 
efficiently treated using the diffusion approximation. How- 
ever many astrophysical problems involve transport through 
substantial regions of optically thin material, and here we 
conduct a test of the MC algorithm that incorporates both 
optically thin and opaque material. We also include a time 
varying source of photons which allows us both to test how 
well the algoritfim retains the coherence of the radiation 
field but also how well it captures the heating and cooling 
of material. 



We adopt a f cm box with a density of 2 x 10 ^ g cm 



and 
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Figure 3. Heat kernel test for the pure scattering case. In the 
upper panel the solution to the diffusion equation is plotted (solid 
lines) for times* = 10"", 2x 10"", lO-^'', and 2xl0~^0s, while 
the radiation energy densities for the same times found using the 
MC algorithm are plotted as crosses. The lower panel shows the 
absolute differences between the analytical and MC-based energy 
densities at the same timesteps (solid, dotted, dashed, and dot- 
dashed lines respectively). 



Figure 4. Heat kernel test for the pure absorption case. In the 
upper panel the solution to the diffusion equation is plotted (solid 
lines) for times t = 10"", 2x 10"", 10'^°, and 2x10-1'' s, while 
the radiation energy densities for the same times found using the 
MC algorithm are plotted as crosses. The lower panel shows the 
absolute differences between the analytical and MC-based energy 
densities at the same timesteps (solid, dotted, dashed, and dot- 
dashed lines respectively). 



10" 



if a; < 0.5 
if X > 0.5 



(27) 



Photons are injected into the left-hand boundary of the 
box from a source of luminosity that varies sinusoidally be- 
tween zero and 10^'^ ergs-^ on a period of IQ-^^ s. Photons 
propogate through the optically-thin part of the box un- 
til they encounter the optically thick-material and are ab- 
sorbed, thus heating the right hand half of the box. This 
heated material emits radiation in an attempt to come into 
thermal balance; since the optical depth is smaller in the —x 
direction these photons are preferentially emitted in that di- 
rection (accompanied by a transport by diffusion in the +x 
direction). Outflow boundary conditions are imposed either 
end of the box. 

Several timesteps from this test calculation are pre- 
sented in Figure [5] Initially the sinusoidally varying radi- 
ation fleld can be seen to propagate at speed c through the 
left-hand half of the box. The optically thick material is 
heated and re-emits radiation which can then be seen prop- 
agating to the left with a sinusoidally varying flux (the am- 
plitude of the variability is lower than that of the impinging 
radiation field, never reaching zero, due to thermal capacity 
of the optically thick material). The diffusion of radiation 
through the optically thick material is also apparent. 

We have checked the energy conservation of the algo- 
rithm as the source radiation is transported, absorbed and 
remitted. In Figure |6] shows the integrated thermal and ra- 
diation energy presented in the box at each timestep (prior 
to any radiation escaping the box at the boundaries). The 
total energy is conserved to better than 1% at all timesteps. 



3.4 Implementation 

The above 1-d, grey tests indicate that the algorithm is fun- 
damentally sound and can reproduce analytical results. We 
therefore implemented version of the algorithm as a mod- 



ule within the TORUS radiative transfer code ( Harries|2000 
Kurosawa et al.||2004| |Pinte et al.||2009| ). The code is writ 



ten in Fortran 90 and follows the radiative transfer on an 
adaptive mesh stored as a quad-tree (2-d) or oct-tree (3-d). 

Since each Monte-Carlo photon packet is essentially an 
independent event the code is straightforwardly parallelized 
under the Message Passing Interface (MPI). Each thread 
holds a copy of the grid in memory, and the work of the 
photon loop is distributed over all the MPI threads. At the 
end of the photon packet loop the sums in equations |10| 
and [11] are made over all threads, and thus estimates for 
the photon energy density and absorption rate are found for 
each cell in the grid, and finally the results are distributed 
back to all threads. The communication overhead for the 
last step is minimal compared to the calculations for the 
photon loop, and the calculation is thus CPU rather than 
thread-communication limited. 

Example CPU times are susceptible to rapid erosion in 
usefulness due to Moore's law, but for the benefit of con- 
temporary readers the science grade calculation described 
in section |4] was run on 16 theads on the University of Ex- 
eter's SGI Altix ICE 8200 supercomputer and each of the 
later timesteps (when the photon stack size had maximized) 
took approximately 4 minutes, of which 20 seconds were 
inter-thread communication. It was found that the calcula- 
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Figure 5. An optically thin/thick radiation propogation test with 
a varying source. The panels (top to bottom) are snapshots at 
t = 10-", 2 X 10-", 4 X 10"" and 6 X 10-" s. Plotted are 
the radiation energy density of the source photons (solid line), 
the thermal energy density of the material (dashed line), and 
the radiation energy density of photons emitted by the material 
(dotted line). The expected maximum extent of the radiation field 
from the source at t = lO-^'^ s is plotted as a dot-dashed line in 
the upper panel. 



1.8x10' 
l.Sxlo' 



1.2x10' 

U IxlO' 




Time (seconds) 



Figure 6. Total energies for the the optically thin/thick radia- 
tion propogation test with a varying source. In the upper panel 
the expected total energy in the box is plotted as a solid line, 
along with the total energy calculated from the MC algorithm 
(dotted line, very close to the solid line), the energy in the radi- 
ation field (dot-dot-dash line), the thermal energy (dashed line), 
and the radiation density emitted from the heated material (dot- 
dash line). The vertical dashed line indicates the expected time 
that the source radiation field enters the optically thick material 
(1.67 X 10^^^ s). The lower panel shows the fraction difference 
between the expected total energy in the box and that derived 
from the MC-mcthod. 



tion scaled almost perfectly up to 64 threads (the maximum 
that we tested). 



benchmark midplane temperature distribution is found after 
~10i«s. 



3.5 Benchmark protostellar disc 

As an initial test we used the 2-d RT benchmark disc de- 
scribed by Pascucci et al. ( 2004[ ). By following the radiation- 
transport over a sufficiently long duration the disc (initially 
at T = OK) will come into radiative equilibrium with the 
central star. This represents a rigorous test of the new al- 
gorithm since the circumstellar material contains both op- 
tically thin and optically thick regions and thus simultane- 
ously incorporates both free-streaming radiation and diffus- 
ing photons. Furthermore the polychromatic nature of the 
algorithm is tested with excess luminosity from the pho- 
tosphere at short wavelengths heating the disc and being 
re-emitted in the mid-IR and sub-mm regimes. 

We adopted an initial timestep of 100 s, and used 10^ 
photons per timestep. In Figure[7]the gradual heating of the 
disc is shown. The rarified regions above and below the disc 
quickly reach thermal equilibrium, while the disc midplane, 
which has a characteristic radial optical depth of ~ ICQ at 
5500 A, approaches equilibrium much more slowly. After ~ 
3 X 10^" s the entire disc has reached a steady state. 

A quantitative comparison with the benchmark solution 
is shown in Figure [S] Good agreement with the published 



4 APPLICATION: PROTOSTELLAR DISC 
VARIABILITY 

The accretion rate onto protostars is inherently variable on a 
wide variety of timescales from hours to years (e.g. [Bouvier 
et al.|2007 Nguyen et al.|2009 1 . This variable accretion flux 
(emitted primarly at wavelengths less than 4000A) is scat- 
tered and reprocessed in the disc to be emitted in the near- 
and mid-IR. There will be an intrinsic delay between the 
optical variability and the disc's IR response, in part due to 
the light-travel time from the star to the disc, but also due 
to thermal lag within the disc itself. Since the temperature 
in the disc decreases radially, one expects the time lag (with 
repect to the optical flux variability) to increase with wave- 
length. Hence by obtaining simultaneous timeseries photom- 
etry in the optical and IR one may attempt to map the 
disc emission using the wavelength-dependent lag, in much 
the same way as reverberation mapping is used to map the 
broad-line region around AGN from the emission line/UV 
continuum variability (e.g. Denney et al. 20091. Of course 
there are a host of other factors that might complicate this 
rather naive interpretation of the variability, such as changes 
to the disc structure, and rotational modulation of either an 
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Figure 7. The Pascucci benchmark disc heated from T = K to thermal equihbrium using the MC algorithm. The temperature of the 
disc is plotted as a logarithmic colour-scale scaled between 1 and 100 K. Distances are in AU. Four snapshots are shown: t = 3.1 X 10* s 
(top left panel), t = 1.02 X lO** s (top right), t = 3.3 X 10^ s (bottom left), and t = 1.05 X lO'' s (bottom right). 



azimuthally structured inner edge (warps) or an anisotropic 
radiation field (hot spots), but it is nonetiieless interesting 
to examine the expected timescales from a simple thermal 
response model. 



4.1 Model parameters 

We adopt a central protostar of radius R* = 2R0 and a 
blacicbody photospheric flux distribution at T^s = 4000 K. 



The recent study by Nguyen et al. ( 2009 1 showed that the 



accretion rate onto a typical protostar varies by a factor 
of ~ 2 on a timescale of days to weeks. We assume the 
star is accreting at a rate that varies sinusoidally between 
5 X lO'^Mgyr"^ and 1 x lO'^Mgyr"^ with a period of 
1 h. This period is comparable with the flushing timescale 
of the magnetosphere and line profile variability on such 



timescales has been observed e.g. Smith et al. ( 1999 1; Kuro- 



sawa et al. (20051. Such rapid variability corresponds to a 



change on a much shorter timescale than the canonical rota- 
tion period, making it easier observationally to distinguish 
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Distance (au) 

Figure 8. The disc midplane temperature for the Pascucci bench- 
mark disc. The heavy line shows the benchmark temperature dis- 
tribution, while the other lines show (from the left) snapshots at 
t = 3.1 X 10" s, t = 1.02 X 10« s, t = 3.3x lO"^ s, t = 1.05 x 10^ s 
and t = 3.3 X lO"' s. 



disc reprocessing effects from rotationally-modulated disc- 
structure effects. 

In order to mimic the additional fiux associated witli the 
accretion we add tfie accretion luminosity as a blackblody 
with a characteristic temperature found by assuming that 
the accretion power is emitted from an area equivalent to 
5% of the stellar photosphere. 

A simple flared structure is adopted for the disc, viz: 



p{r,z) 



■ por exp 



2 h{ry 



where 

h = ho{r/r())^ 



(28) 



(29) 



with a = 2.25 and /3 = 1.25. The disc scale height is set to 
ho — 125 au at rg = 100 au. In order to simulate the trunca- 
tion of the inner disc by the magnetosphere we assume an 
inner hole of 10 R,, while the outer radius is 300 au. We fix 
Po by assuming a disc mass of 0.01 Mq. The disc is assumed 
to contain the dust size distribution and chemisty described 



by Kurosawa et al. (20041, and for the purposes of this test 



we assume isotropic scattering. 



4.2 Method 

The disc is initially brought into radiative equilibrium with 
the central object (assuming a constant accretion rate of 5 x 
10~^ M0 yr~^) by using a time-independent algorithm. The 
time-dependent method is then turned on (which defines 
t = s) and the RT is followed for ten periods at a timestep 
of 18 seconds. This short timestep allows us to adequately 
resolve the shortest timescales that are of interest, which is 
basically the light-crossing time from the central object to 
the inner disc (47 s). 

Since the photon packet estimate of the radiation field 
is only defined within a sphere of radius ct around the cen- 
tral object it is only cells within this radius whose thermal 
properties are changed at each timestep, with the rest of 
the disc held in thermal equilibrium. Free-streaming photon 



packets that have a negligble probability of interacting with 
the disc (found by integrating the optical depth along the 
packet's path) are deleted from the stack at the end of each 
timestep. This has the advantage of significantly reducing 
the memory requirement of the photon stack. 

The shape of the cross correlation function should be a 
function of the system's viewing angle (inclination) as dif- 
ferent parts of the disc have different star-disc-observer path 
lengths. In the following we adopt an inclination of 60° when 
calculating the SEDs. 

Initially a set of photon packets were generated from the 
volume of grid outside the ct sphere defined by the length 
of the timeseries. Since this volume is at constant temper- 
ature it was possible to calculate this contribution to the 
SED in a time-independent manner. Subsequently at each 
timestep 10^ new photon packets (photospheric or thermal 
disc) were generated, and these photons were tagged by their 
generation time. A 'peel-off' technique was then employed: 
The light travel time and optical depths to the observer 
were calculated, along with the probability that the pho- 
ton packet was emitted towards the observer direction. The 
photon packets were then binned in the observer's frame, 
both spectrally and temporally. The new photon packets, 
plus the photon packet stack, where then followed for a sin- 
gle timestep, with packet peel-off and binning occurring at 
each scattering event. 

4.3 Results 

We have computed the cross-correlation function (CCF) of 
monochromatic time-series from 3000 A to 10 fim against a 
fiducial time-series at 3000A (see Figure l9| . It is immedi- 
ately apparent that the lag increases with wavelength due 
to a combination of hght travel time to the disc, and the 
thermal lag as the increased radiation field heats the lo- 
cal disc material. The maximum value of the cross correla- 
tion function decreases at longer wavelengths as the thermal 
timescale of the disc begins to dominate over the variability 
timescale. 

I have quantified the lag by fitting a Gaussian to each 
CCF peak as a function of wavelength. The central positions 
of these Gaussians are plotted in Figure [To) The lags at blue 
wavelengths (^ 2 /im) are short, in fact less then the hght 
travel time between the star and the disc, this is principally 
because there is significant direct photospheric emission at 
these wavelengths (with zero lag) in addition to the scattered 
and reprocessed components. Redwards of ~ 2 fim the flux 
is dominated by thermal emission from the disc, and we 
see that the lag increases monotonically (within the errors 
on the CCF Gaussian fits) with wavelength (although the 
strength of the correlation is decreasing). 

Observing this correlated phenomenon would be a chal- 
lenging proposition, necessitating high cadence temporal 
sampling in the optical and near-IR simultaneously. Opti- 
mally one would select ground-based UBV photometry com- 
bined with Spitzer IRAC photometry. Although it may be 
possible to conduct the IR observations into the mid-IR it is 
clear that the correlation is dropping rapidly at these wave- 
lengths, and it is not at all clear the necessary temporal 
sampling could be achieved using ground-based facilities. It 
may be that systems with a larger inner-hole, such as Herbig 
AeBe stars ( Monnier et al.|2005 l would be superior targets 



Time- dependent radiative transfer 9 



Time lag (seconds) 



Figure 9. A cross-correlation image of the monochromatic disc 
hghtcurve at 3000 A against wavelengths from 3000 A to 10 fira. 
The image is a linear greyscale of cross-correlation function scaled 
between 0.8 (white) and 1 (black). The bottom row of the image 
is the autocorrelation function of the 3000 A lightcurve. 
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of sources embedded in an arbitrary distribution of ab- 
sorbing, emitting and scattering material. The new algo- 
rithm, based on the MC radiative-equilibrium method of 



Lucy ( 1999 1 can be used in 3-dimensions and scales almost 



perfectly under parallelization. It has advantages over the 
FLD approximation in that it is polychromatic and correctly 
treats the directionality and flux of the radiation field in the 
optically thin limit. Note that although the algorithm as pre- 
sented here is only applicable under conditions of LTE, the 
method may be straightforwardly extended to the non-LTE 
regime e.g. 



Carciofi & Bjorkman ( 2006 1 



I have applied the new method to the problem of a cir- 
cumstellar disc illuminated by a protostar with a periodic 
time-variable accretion rate. I have shown that the lag be- 
tween the blue continuum resulting from the accretion hot 
spots and the reprocessed IR radiation from the disc is a 
strong function of wavelength. It appears that photometric 
time-series data in the blue part of the optical spectrum, 
combined with equally intensive near-IR time-series could 
be used to probe the geometrical and thermal structure of 
the disc, although it is likely that complications arising from 
stochatic variations in the disc structure close to the inner 
edge could mask the correlation calculated here. 

In future the time-dependent method will be incor- 
porated in radiation-hydrodynamics calculations of proto- 
stellar disc fragmentation, extending the time-independent 



transfer simulations of Acreman et al. (20101. It is well es 



tablished that the thermal properties of the disc strongly 



affect the fragmentation (Boss 2008 Stamatellos & Whit- 
worth|[2008 1 , and it is also clear that the FLD approxima- 
tion does not capture the full physics of the transport in the 
disc itself. 
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Figure 10. The time lag with respect to the 3000A lightcurve 
plotted as a function of wavelength (filled circles) . The light travel 
time between the stellar centre and the inner disc edge is plotted 
as a dashed-line. 

to Classical T Tauri stars, allowing light-travel time effects 
to remain dominant at longer timescales. 

Of course in reality the underlying accretion luminosity 
is unlikely to vary periodically, but more stochastic varia- 
tions could still be investigated. It may be possible to dis- 
tinguish between disc thermal responses and scattering by 
observing polarized light, which will be dominated by the 
scattered light contribution. 



5 CONCLUSIONS 

I have presented a new method for computing time- 
dependent radiation transport for an arbitrary distribution 
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